Virtual greenhouse model

Niels Holst, Aarhus University, Denmark, 25 May 2023.

R scripts

The documentation comes with R script that demonstrate the equations and logic of the model. Before you run any of the scripts, you need to change the R work folder. This is defined in the setup.R script, which begins like this:

Model structure

The model is updated iteratively in a loop at a time step of Δt (s), typically in a range of 2 to 5 minutes; however, the energy and water budgets may be resolved at a finer time step as needed (dividing Δt into n sub-steps). The model runs from the chosen start date until the end date; however, beneath the surface it will always start 24 hours before the start date to allow model state variables to evolve from their initial values and reach realistic values.

Model components are updated in the following order:

  1. Calender updates current date and time.

  2. Outdoors updates climatic variables based on the weather file. Readings from the weather file are interpolated linearly to estimate values for the current date and time. If the weather file does not contain a Tsky column, sky temperature is estimated by a model. Soil temperature is currently roughly modelled (this needs an update).

  3. Greenhouse construction contains the dimensions of the greenhouse and a description of each of the six greenhouse faces (traditionally, glass, see cover layer).

  4. Setpoints are defined for the green house actuators, such as growth lights, heat pipes, ventilation, screens and CO2 diffusers, to keep indoors climate within the desired bounds.

  5. Controllers use the setpoints to compute the settings of actuators.

  6. Actuators represent growth lights, heat pipes, ventilation, screens, CO2 diffusers, etc., which control the greenhouse climate. They are governed by the setpoints via the controllers.

  7. Indoors holds the greenhouse compartments of layers and volumes and keeps their state variables updated, reflecting the current outdoors driving variables and actuator settings.

  8. Output writes the chosen variables to the output file.

Boxscript example with calendar

The simulation begins 24 hours before the declared beginning of the simulation. This allows model state variables to attain sensible initial values. Here is a possible layout of a boxscript. The lonely : characters denote left-out lines:

A global box is declared to hold global settings, i.e., parameters that serve as inputs to one or more of the boxes further down in the boxscript. A beginDate port is created (the ampersand & means that a new port is created) to start the simulation on 1 August 2022. The calendar actually starts one day before, as set by its begin port (line 9); the initial period . denotes an existing input port. At the very end of the boxscript, an OutputSelector box asks for the output to be written starting on 1 August (line 35). In effect, the initial extra day's output will be discarded. The simulation will run until the final date end set by the calendar box; the time step has been set to 15 minutes (lines 10-11).

Greenhouse virtualization

Greenhouse compartments

Greenhouse dimensions

The greenhouse consists of n spans each of width w (m), length (m) and wall height h (m). Each span has a roof on top angled at u (°).

From these parameters we can derive several measures of greenhouse geometry. The area of a single face, either side, end or roof, is

(1)Aside=h(2)Aend=wh+w24sin(u)(3)Aroof=w21+sin2(u)

from which follows the total area of side, end and roof faces, and the total cover area (all in m2),

(4)Aside=2Aside(5)Aend=2nAend(6)Aroof=2nAroof(7)Acover=Aside+Aend+Aroof

The ground area (m2) enclosed by the greenhouse is

Agh=nw

and the volume contained (m3) is

(8)Vgh=nAend

Further useful measures are the cover/area ratio Cai (m2 cover/m2 ground),

(9)Cai=AcoverAgh

and the volumetric area index Vai (m) (also known as the average height),

(10)Vai=VghAgh

Furthermore, the density of air ρair=1.19 kg/m3 and it heat capacity Cair=1020 J/kg/K leads to the heat capacity of indoors air Cgh (J/K/m2 ground),

(11)Cgh=ρairCairVai

The equations above have been coded in the greenhouse-geometry.R script::

Boxscript example

In boxscripts, greenhouse dimensions are set in the geometry box inside construction:

Volumes

The greenhouse model holds three volumes: indoors, outdoors and soil. The indoors holds state variables, while the outdoors and soil hold driving variables. All three hold temperature T (°C), while indoors and outdoors hold relative humidity H% (%) and CO2 concentration [CO2]ppm as well. In some relations other units are more convenient, e.g. absolute humidity Habs (kg/m3), water vapour pressure Hpres (Pa), or absolute CO2 concentration [CO2]abs (g/m3). The variable are subscripted to denote the volume, e.g., for temperature Tin, Tout and Tsoil.

Volumes may exchange heat and humidity with neighbouring volumes and layers.

Layers

Besides layers the greenhouse and its environment is conceptualised as a stack of i=1..n horizontal layers. This is a typical stack:

  1. sky

  2. cover (e.g., glass)

  3. screen

  4. growth light

  5. plant canopy

  6. heat pipe

  7. floor

Depending on the scenario, layers can be removed or added (e.g., additional screen layers or a table layer ). All layers, whether they are horisontal or inclined, uniform or heterogeneous, thick or thin, are simplified to describe an average 1m2 ground area of the greenhouse. For screens this area will vary according to what degree they are drawn. Where the unit m2 is used, it always refers to m2 ground. Sometimes it is spelled out "m2 ground" to avoid confusion.

The radiative properties of each layer is parameterised for three intervals of the spectrum:

PAR is a subset of the short-wave radiation. The PAR budget is resolved only to compute how much is absorbed by the plant canopy. Its quantum units reflect this purpose. In contrast, the short- and long-wave radiation both contribute to the radiative exchange of heat among layers. Remember that W = J/s.

energy-budget-layers

 

Each layer i is characterised by its radiative properties: absorptivity (αir), reflectivity (ρir) and transmissivity (τir) for radiation coming from above and similarly for radiation coming from below (αir,ρir,τir), all for each type of radiation r (PAR, sw, lw). The radiative parameters always add up to unity, αir+ρir+τir=1 and αir+ρir+τir=1. We apply the common engineering assumption that emissivity equals the absorptivity for long-wave radiation, i.e. for all layers ϵi=αilw and ϵi=αilw.

The superscript r was left out of the figure above for simplification but it is important, as it allows different parameter settings for different radiation intervals. This is essential, e.g. for the distinct difference in transmissivity of glass for short- vs. long-wave radiation. Moreover, the separate parameter sets for down- and upgoing radiation allow for asymmetry, e.g., of screens.

Layers are further described by their area-specific heat capacity (Ci; J/K/m2 ground) and temperature (Ti; °C).

Boxscript example with cover and screens

The construction box holds two child boxes, geometry and shelter. The shelter box holds definitions of different cover and screen materials inside the covers and screens box, respectively. In the boxscript below, one kind of cover is defined (named glass) and four kinds of screens (named energy, shade1, shade2 and blackout). You can name the materials as you want.

The faces box (lines 67-98) holds the named materials of cover and screens for each of the six shelter faces. There may be more than one layer of screens installed. These screen layers are written in sequence starting nearest the glass. For example, with up to three layers installed on a face, you might have faces with screens specified as energy+shade1+blackout, none+shade2+blackout, energy+none+blackout or none+none+blackout, or any other combination of screen materials defined by the names of the shelter/screens/* boxes. All six faces must have the same number of screens defined. If a face has no screen in a certain layer, it is written as none.

Any layer has 11 input ports defining its basic radiative properties (lines 14-21), U-values (lines 22-23) and heat capacity (line 24). Cover boxes holds two additional inputs (lines 25-26) specifying the current status of whitening ('chalk'), which is provided by the chalk actuator. In a similar vein, Screen boxes hold one additional input state, which refers to a screen actuator. From the code above, you can infer that in this model there are three screen actuators:

In addition to cover and screens, each face has a weight parameter that is used when calculating the average radiative properties of the greenhouse cover and each screen layer. The settings above suggest that side2 is north-facing.

Processes of heat exchange

Radiative exchange

Radiation might be emitted by a layer, both down- and upwards, in terms of PAR (Eipar,Eipar; μmol/m2/s) and short- and long-wave radiation (Eisw,Eisw,Eilw,Eilw; W/m2). For some layers (sky, cover, screen, floor) , the emission of long-wave radiation is calculated from their temperature, according to the Stefan-Boltzmann Law,

(12)Eilw=σϵi(Ti+T0)4Eilw=σϵi(Ti+T0)4

where T0=273.15 °C converts to absolute temperature. For other layers (growth light, plant canopy, heat pipe), the emission is calculated by specific sub-models..

The emitted radiation is intercepted by neighbouring layers, its ultimate fate determined by the characteristics (absorptivity, reflectivity, transmissivity) of the whole stack of layers, as worked out in the energy budget. The net radiative absorption of a layer (i.e. its absorption minus emission) is denoted by (Aipar,Aipar; μmol/m2/s) and (Aisw,Aisw,Ailw,Ailw; W/m2).

Convection and conduction

The exchange of heat between layer i's upper and lower surface and the confronting volume is described by Ui on the upper side and Ui on the lower side (both in W/K/m2). The mechanism of heat exchange is convection for air volumes and conduction for the soil volume but despite that they are all superscripted as Aiconv.

The convective/conductive heat absorbed by layer i (W/m2) from the volume above and below is, respectively,

(13)Aiconv=Ui(TViTi)Aiconv=Ui(TViTi)

where the temperature of the volume above and below layer i is denoted by TVi and TVi, respectively.

Advection

Ventilation rate integration

Advective heat transfer occurs when outdoors air replaces indoors air. This is due to intentional ventilation as well as inadvertent leakage. We formulate the ventilation rate as, "how many times greenhouse air is replaced per hour" v (h1).

During a time step Δt the indoors temperature Tin and absolute humidity Hinabs will approach the outdoors Tout and Houtabs. Considering either variable and assuming a perfect mixture of air inside the greenhouse, the rate of change in the indoors variable yin is

dyindt=v(youtyin)

Integration gives us

(14)yin(t)=yout[youtyin(0)]exp(vt)

The relation is illustrated by the advection.R script with Tout at 20°C and Tin(0) either 15°C or 25°C. The ventilation rate was set at v=3 h1.

To calculate the change in indoors temperature ΔTin during a time step Δt (s) eq. 14 can be rewritten as

(15)ΔTin=[ToutTin(0)][1exp(v3600 s/hΔt)]

With a time step of, for example Δt=2 minutes, the decrease (blue curve) or increase (red curve) in indoors temperature would be ΔTin=0.48 K (by eq. 15). Assuming a volumetric area index of Vai=3.94 m3/m2 (eq. 10) and with air heat capacity Cair=1214 J/m3/K, we get the following energy flux (positive in the case of the red curve),

(16)Ainadv=0.48 K1214Jm3K1120 s3.94 m=19.0Wm2

The change in indoors absolute humidity ΔHinabs can be arrived at in a similar fashion.

Eq. 15 is an example of a local integral; it is solved assuming everything else remains constant. It is applied here because high values of v could make a straightforward Euler integration imprecise. In this case, the difference is not so great. Euler integration would have yielded

ΔTin=[ToutTin(0)]v3600 s/hΔt=0.50 K
Maximum ventilation rate

The maximum ventilation rate vmax (h1) is obtained when the vents are fully open. What this ventilation rate amounts to depends on the wind speed w (m/s) and the difference between the indoors Tin and outdoors Tout temperature (°C):

(17)vmax=AventAgh(cww)2+(cTTinTout0)2

with discharge coefficients for wind cw (s/m/h) and temperature cT (/K/h). The total discharge is found by vector addition, considering the two forces operating independently (orthogonally). It is applied to the relative vent area; Avent (m2) is the total area of vents and Agh (m2) is the greenhouse floor area.

Latent heat

Latent heat is related to the phase change between water in the gas and the liquid phase. The phase change costs (by evaporation) or releases (by condensation) energy according to the latent heat of water vapourisation λ=2454 kJ/kg.

In the model, transpiration cools the canopy and condensation heats the cover. The transpired water adds to air humidity, which is then considered to have increased its content of latent heat. Water vapour bears latent heat in the sense that if it condenses, the heat used for the initial phase change from liquid to gas, will be used in the phase change back from gas to liquid.

The model does not include the cooling effect of water evaporating from wet surfaces, nor the wetting and drying of screens.

Layer characterization

Sky layer

The radiation emitted downwards from the sky includes Eskypar, Eskysw and Eskylw. PAR and short wave radiation originates from the sun, while long wave radiation is emitted from the sky as such.

While Eskysw is measured routinely, this is not the case for Eskypar. It can be approximated by Eskypar=0.5Eskysw.

Sky emissivity ϵsky (and thereby sky absorptivity αsky) is estimated from the dew point temperature Tdew (°C) by the empirical relation (ref.)

ϵsky=0.732+0.00635Tdew1

Sky temperature Tsky (°C) is then derived from outdoors temperature Tout (°C) according to Karn et al. 2019

Tsky=ϵsky4(Tout+T0)T0

where T0=273.15 °C.

Sky reflectivity ρsky is zero and consequently transmissivity is τsky=1ϵsky.

Cover layer

The cover layer is commonly made out of glass with the following typical characteristics:

Glass propertiesPARShort-waveLong-wave
α,α,ϵ,ϵ0.10.10.9
ρ,ρ0.10.10.1
τ,τ0.80.80.0

In general, the six faces of the greenhouse may carry different covers with different radiative properties. For the cover layer as a whole, each radiative parameter is summarised as the average over all faces, weighed by face area.

The density of window glass is typically 2.5 kg/m2/mm, which with a thickness of 4 mm gives a density of 10 kg/m2 cover. The heat capacity of glass is 840 J/kg/K (Engineering Toolbox, 2023). If we furthermore assume a cover/ground ratio of 1.24 m2 cover/m2 ground (eq. 9), we get a heat capacity of the cover of

(18)Ccover=10 kg/m2 cover1.241 m2 cover/m2 ground840 J/kg/K=10421 J/K/m2 ground

The cover emits long-wave radiation symmetrically down and up according to the Stefan-Boltzmann Law (eq. 12).

The U-values over both the inner and outer surface can be expected to increase with the wind; immobile air is a very good insulator. For the outer surface of the cover, we use the empirical relation,

(19)Uwind=2.8+1.2w0.8

where w is the outdoors wind speed (m/s). The U-value of the inner surface is set to 1.2 W/K/m2.

Cover averaging

Since the six cover faces may differ in their characteristics, we may need to average their parameter values. Some examples applying the equations below of are shown in the cover-averaging.R script.

Cover averaging of radiative parameters

Let's assume that the two end faces differ from the other four faces; we set α=0.15, ρ=0.35 and τ=0.60 at the ends walls, while elsewhere, α=0.10, ρ=0.10 and τ=0.80. In the full model, these values have been corrected for the effects of chalking before the averaging. Furthermore, let's weigh the end walls by wend=0.25 and one side by wside1=0.20 while the rest have full weight wother=1.0.

For the radiative parameters (α, ρ and τ), we get their weighted averages (αcover, ρcover and τcover) for the cover as whole by

(20)α0=i=16wiAiαii=16wiAi,ρ0=i=16wiAiρii=16wiAi,τ0=i=16wiAiτii=16wiAiαcover=α0α0+ρ0+τ0,ρcover=ρ0α0+ρ0+τ0,τcover=τ0α0+ρ0+τ0

where subscript i stands for the six cover faces, and Ai (m2) is the area of each face. Since α0+ρ0+τ0=1 is not guaranteed, these values are subsequently rescaled to a sum of 1.

With the numbers above and the illustrative greenhouse dimensions we have (with parameters pair-wise in the order: roofs, sides, ends)

w¯=(1,1,1,0.20,0.10,0.10)A¯=(5459,5459,350,350,394,394) m2 coverα¯=(0.10,0.10,0.10,0.10,0.15,0.15)ρ¯=(0.10,0.10,0.10,0.10,0.35,0.35)τ¯=(0.80,0.80,0.80,0.80,0.60,0.60)

which gives us the averaged values,

αcover=0.101ρcover=0.104τcover=0.795

The roof faces are dominating by their size; hence the average parameter values are close to the roof values. If we change all weights to one wi=1, we get a larger impact of the ends: αcover=0.102, ρcover=0.115 and αcover=0.783.

Cover averaging of heat capacity

Heat capacity is averaged to apply to the greenhouse ground area (eq. 18). If, for example, the largest area is covered by glass with a heat capacity of Cglass=8400 J/K/m2 cover, while the ends have a more sturdy material with Csturdy=18000 J/K/m2 cover, we can calculate the total heat capacity of the cover per greenhouse area Ccover (J/K/m2 ground) as

(21)Ccover=i=16AiCiAgh

where Agh (m2) is the ground area of the greenhouse. We apply no weights to the different greenhouse faces, other than their areas, when averaging the heat capacity.

Using the the illustrative greenhouse dimensions we have

A¯=(5459,5459,350,350,394,394) m2 coverC¯=(8400,8400,8400,8400,18000,18000) J/K/m2 coverAgh=10000 m2 ground

which gives us

Ccover=11177 J/K/m2 ground

As expected, this heat capacity is larger than the one calculated for and all-glass greenhouse Ccover=10421 J/K/m2 ground in eq. 18.

Cover averaging of U-values

The Ui values (W/K/m2 cover) are conductances (one for each face, i=1..6). They can also be looked upon as resistances by taking their inverse, Ri=1/Ui. If we assume that the six resistances work like electrical resistors in parallel, we get the total resistance of the greenhouse cover Rcover as

1Rcover=i=161Ri

It follows that conductances are simply additive. Thus the conductance of the greenhouse cover as a whole is arrived at by the sum of the faces. Relating this to the ground area of the greenhouse Agh (m2 ground), we get

(22)Ucover=1Aghi=16AiUi

As for heat capacity (eq. 21), we apply no weights to the different greenhouse faces, other than their areas.

For example, If the cover material has Uglass=1.2 W/K/m2, except for the end walls which have Upoly=0.3 W/K/m2, we get with the illustrative greenhouse dimensions,

A¯=(5459,5459,350,350,394,394) m2 coverU¯=(1.2,1.2,1.2,1.2,0.3,0.3) W/K/m2 coverAgh=10000 m2 ground

from which follows

Ucover=1.42 W/K/m2 ground

The reason why the average Ucover is larger than any Ui is that the ground area is smaller than the cover area.

Screen layer

Screens are often constructed with asymmetric radiative properties. Here, for example, are the radiative parameters of a polyester screen with an aluminium surface on the upper side (from Bailey 1981):

Screen propertiesLong-wave radiation
α,ϵ (polyester)0.57
ρ (polyester)0.43
α,ϵ (aluminium)0.07
ρ (aluminium)0.93
τ,τ0.0

The values in the table are for long-wave radiation only. They are most likely different for short-wave radiation but documentation for the physical characteristics of screens are hard to come by; manufacturers seem reluctant to provide the information.

Screens may be more or less drawn . If we denote how much a curtain is drawn by its position p[0;1] then the effective radiation parameters are calculates as pα, pρ and 1pαpρ, and likewise for α, ρ and τ. For the screen layer of the greenhouse as a whole, each radiative parameter is summarised as the average over all the six faces, weighed by face area, as for the cover layer. If screens are kept half-drawn for extended periods, or if faces differ much in screen material, this averaging could turn out highly unrealistic. If the greenhouse is equipped with more than one layer of screens, the screens are stacked as separate, independent layers.

Screens may be constructed of composite materials which complicates the calculation of their heat capacity. For the screen exemplified here, we found quite close values of Cpolyester=1.1 J/K/g and CAl=0.89 J/K/g. Assuming that polyester is the dominant material and a screen density of 80 g/m2 screen, we get a heat capacity of the screen of

(23)Cscreen=80 g/m2 screen1.24 m2 screen/m2 ground1.1 J/K/g=109 J/m2 ground

This calculation implies that all six greenhouses faces are equally covered with the same screen. If that is not the case then Sscreen must be found by calculating the heat capacity of each face separately and adding them up to give Sscreen. Note that the calculated heat capacity of the screen is ~100 times less that of a glass cover (eq. 18).

The screen emits long-wave radiation according to the Stefan-Boltzmann Law (eq. 12). With the example above, the radiation would be asymmetric because downwards and upwards emissivity differ; this screen was designed to keep the heat indoors.

Screen averaging

The averaging of screens in a layer follows the same principles as cover averaging, except that screens also have a position p[0;1]. In a given layer, there may be screens on some layers but not on others; these are given full transparency, τ=τ=1.

Thus for the radiative parameter q (i.e. α, ρ, α or ρ), we get the weighted average qscreen for the screen layer as a whole (cf. eq. 20),

qscreen=i=16piwiAiqii=16piwiAi

with transmissivity as in eq. ???,

τscreen=1αscreenρscreenτscreen=1αscreenρscreen

The heat capacity and U-value of screens is assumed to be in effect only when they are drawn. Hence we get (cf. eqs. 21 and 22) the total heat capacity of the screen layer Cscreen (J/K/m2 ground),

Cscreen=i=16piAiCiAgh

and the total U-value of the screen layer Uscreen (W/K/m2 ground),

Uscreen=1Aghi=16piAiUi

Some examples are shown in the screen-averaging.R script.

Screen layering

Screen layers function just like any other layer (see Layers). Thus there is need for special mechanisms to account for screen layering.

Growth light layer

Lamps are installed to provide growth light. Even so, their production of heat is a significant side effect. Light is traditionally emitted downwards from a position above the canopy. If the installed lamp power (including any ballast or driver needed by the lamps) is Plamp (W/m2) and the PAR efficiency is elamp (μmol PAR/J) then the emission of PAR downwards Elamppar and upwards Elamppar (μmol PAR/m2) is

Elamppar=elampPlampElamppar=0

Lamps may loose efficiency with age, which is reflected in a reduced elamp.

Energy is dissipated from the lamps by three routes: short-waved and long-waved radiation and convection. The proportions are denoted plampsw, plamplw and plampconv, where plampsw+plamplw+plampconv=1. While the short-wave radiation is assumed to be heading only downwards, the long-wave radiation and convective heat are both assumed to be emitted equally downward and upwards. Thus we have

Elampsw=plampswPlampElampsw=0Elamplw=Elamplw=plamplwPlamp/2Elampconv=Elampconv=plampconvPlamp/2

The growth light layer has a 100% transmissivity (τ=τ=1) and is attributed neither a temperature nor a heat capacity.

Plant canopy layer

Stanghellini (1987 ,eq. 2.1) sets up the following energy balance for the plant canopy

Rn=H+LE+M+J[W/m2]

The equation states (in her notation) that the energy of net radiation absorbed by the canopy Rn is spent on

Like Stanghellini we will ignore M and J in the energy budget (they both make quite small contributions to the total budget) and follow her logic which leads to, surprisingly, that H does not need to be taken explicitly into account either.

Radiative properties

Goudriaan (1977) provides equations (his eqs. 2.21 and 2.26) to calculate the reflectivity of a plant canopy with leaf area index Lai and extinction coefficient k. The reflection coefficient ρh depends on the scattering coefficient σ,

ρh=11σ1+1σ

σ is a species-dependent parameter commonly set to σ=0.2 (Kropff and Laar 1993, p.38), which gives ρh=0.0557.

The original eq. 2.26 includes the reflectivity of the surface underneath the canopy (the floor in our case). However, we will resolve the distribution of radiation among the greenhouse layers otherwise. Consequently, we can set the reflectivity of the underlying surface in the original equation to zero. Thus we get the reflectivity of the canopy, which is symmetric for up- and down-going radiation, based on his eq.2.26,

(24)ρplantr=ρplantr=exp(krLai)exp(krLai)ρhexp(krLai)exp(krLai)/ρh

where superscript r denotes par, sw or lw radiation. Commonly found values for kpar and ksw are in the range 0.7 to 0.8, while leaves are impenetrable to long-wave radiation resulting in klw=1.

The absorptivity of the canopy is simply

(25)αplantr=αplantr=1exp(krLai)

and the transmissivity

(26)τplantr=τplantr=1ρplantrαplantr=1ρplantrαplantr

The plant-canopy-layer.R script shows how the parameters change with increasing leaf area index for short-wave (k=0.7) and long-wave radiation (k=1):

Canopy temperature

Note: The figures in this sub-section were generated by the 5-energy-budget-plant.R script.

The model of canopy temperature is taken from Stanghellini (1987 ,eq. 3.5):

(27)Tplant=Tin+ri+re2LaiρaCaAplantrad1γ(HsatpresHinpres)1+δγ+rire

with

The saturated vapour pressure Hsat(T) as a function of temperature is described by the empirical equation (FAO, 2023, eq. 11),

(28)Hsat(T)=610.8exp(17.27TT+237.3)

It increases with temperature, as seen in the figure:

The slope δ(T) as a function of temperature is described by the empirical equation (FAO, 2023, eq. 13),

(29)δ(T)=4098(T+237.3)2Hsat(T)

The slope is increasing with temperature too, as seen in the figure:

The internal leaf (stomatal) resistance against water vapour riH2O is computed by the Ball-Barry (ref.) function,

(30)riH2O=[g0+g1Hinrel100%PncinCO2]1

with species-specific coefficients that we here set to g0=0.1 m/s and g1=1.64 m3/mol. Other variables are

Here is an example, how riH2O behaves with increasing humidity at various CO2 levels and with Pn=2.0 μmol CO2 / m2 leaf / s. There is not much variation from the maximum set by g01:

The external leaf (boundary layer) resistance to water vapour reH2O can be computed according to Stanghellini (1987, eq. 2.53):

(31)reH2O=1174|TplantTin|+207u24

where characterizes leaf dimensions and u (m/s) is air velocity at the leaf surface. The constants 1174 and 207 were estimated empirically in a tomato crop by Stanghellini.

This figure shows the behaviour of reH2O with increasing wind speed, at different leaf temperatures and leaf dimensions with Tin=24°C:

Both and u are difficult to estimate. should be smaller for smaller leaves. u should increase with the greenhouse ventilation rate (or if rotors are running) and decrease with Lai, as a denser canopy will reduce the average wind speed at leaf surfaces.

Note: Due to all the uncertainties, we use a fixed value in the model defaulting to

(32)reH2O=200 s/m

Everything above taken together, we get the following linear response of plant temperature Tplant (eq. 27) to the radiation absorbed by canopy Aplantrad with Lai=1.9. The indoors temperature was set at Tin=25°C and is shown as a black line in the figure:

Clearly, a high air humidity works against plant transpiration and leads to a higher leaf temperature, even above the ambient temperature of the greenhouse: Relatively cool leaves are a sign of a sound microclimate.

Transpiration

Note: The figures in this sub-section were generated by the 5-energy-budget-plant.R script.

The model of canopy transpiration rate Hplanttrans (kg/m2 ground/s) is taken from Stanghellini (1987 ,eq. 3.4):

(33)Hplanttrans=δγAplantrad+2LaiρaCaγre(HsatpresHinpres)λ(1+δγ+rire)

with symbols defined as for eq. 27. With the same parameter settings as for the previous figure, we get this transpiration rate of the canopy (shown in units of g/m2 ground/min):

Obviously, transpiration is larger in a drier climate and it increases with the absorbed radiation. Note that leaf temperature (previous figure) too increased with absorbed radiation. This means that the increase in transpiration with increased absorbed radiation was not sufficient to keep leaf temperature constant; that would have acquired an even higher transpiration.

The phase change from water (in the plant) to water vapour (in the greenhouse air) causes an increase in the latent heat carried by the air. This heat, which is equal to λHplanttrans, is shown in the figure below,

The 1:1 dashed line has been added to highlight the balance between the radiation absorbed by the canopy and the latent heat lost from the canopy. Below the line the canopy is warmer than the air, above the line the canopy is cooler than the air. Compare this with the figure of canopy temperature above. At a relative humidity of 90%, both the blue lines cross the respective reference lines at ca. 110 W/m2. This verifies the consistency of eqs. 27 and 33. The radiative energy absorbed by the canopy ends up in latent heat bound in the air and in the heating of the canopy itself.

The convective model of heat transfer does not match the complex physiology and physics of the canopy layer; all heat fluxes are taken care of by eqs. 27 and 33. Therefore the canopy layer is given a U-value of zero to make it fit into the stack of generic layers.

Heat pipe layer

Heat pipes are installed to heat the greenhouse by convection and long-wave radiation. The transmissivity of the heat pipe layer is τpipe=τpipe=1 for all wave lengths.

The drop in temperature (ΔTpipe; K) from the inlet to the outlet is modelled by the empirical equation

(34)ΔTpipe=(TpipeinletTin)[k(b1)Δtpipe+(TpipeinletTin)1b]11bwhenTpipeinlet>Tin

where k and b are parameters calibrated to the greenhouse, Δtpipe (s) is the transit time of pipe water, Tpipeinlet (°C) is the water temperature at the pipe inlet, and Tin (°C) is the greenhouse air temperature. If the inlet temperature is not above the greenhouse temperature then ΔTpipe=0.

The transit time Δtpipe (min) is related to the flow rate v˙pipe (m3/h) and pipe volume Vpipe (m3) as

Δtpipe=Vpipev˙pipe

As an example, a pipe with an inner diameter of 30 mm installed at a density of 2 m/m2 in a greenhouse area of 10 000 m2 holds the volume,

(35)Vpipe=π4(30 mm1 m1000 mm)22mm210000 m2=14.14 m3

which at a flow rate of v˙pipe=20 m3/h gives a transit time of

Δtpipe=14.14 m320 m3/h60 minh=42.4 min

The energy lost from the heat pipe Epipetot (W/m2 ground) is related to the temperature drop through the pipe ΔTpipe and the flow rate v˙pipe,

(36)Epipetot=CwaterΔTpipev˙pipeAgh1000 kgm31 h3600 s

where Cwater=4184 J/K/kg is the heat capacity of water and Agh (m2 ground) is the ground area of the greenhouse. Proper conversion of units must be observed. With ΔTpipe=30 K, for example, we get

Epipetot=4184 J/K/kg30 K20 m3/h10000 m21000 kgm31 h3600 s=69.7 W/m2

As an example, we set k=0.0063 and b=1.25 with an indoors temperature Tin=20°C. We choose a range of inlet temperatures Tpipeinlet=20..80°C and two different flow rates v˙=(10,20) m3/h. Other parameters keep the values used in the examples above. We then get the following drop in water temperature ΔTpipe (eq. 34) reached at the pipe outlet and the associated energy Epipetot (eq. 36) lost to the greenhouse (plots generated by the heat-pipes.R script):

Heat pipe energy is emitted as a combination of radiant and convective heat. We denote the proportion of long-wave radiation ppipelw[0;1] and, assuming that it is radiated equally down- and upwards, we get

Epipelw=Epipelw=ppipelwEpipetot/2

The other part is conveyed to the greenhouse air. Since convective heat is defined relative to the layers (eq. 13), the convective heat flow from the heat pipes is negative,

Apipeheat=(1ppipelw)Epipetot

Since the volume above and below the heat pipes is the same, i.e. the greenhouse air, we let Apipeheat account for the whole heat exchange and set Apipeheat=0 (see eq. 13).

Floor layer

A floor made of concrete has typical values αfloor=0.6 and ρfloor=0.4. It emits long-wave radiation upwards according to the Stefan-Boltzmann Law (eq. 12) and exchanges heat with the soil by conduction (eq. 13). We assume a good insulation of the floor against the soil setting Ufloor=0.1 W/K/m2.

Soil temperature

Soil is not treated as a layer in the model but as a volume below the floor layer (see Volumes and Layers). Anyway, the model of soil temperature is naturally presented here.

Soil is defined only vaguely by the model. As a volume it is endless in reality, and it has got not a single temperature but a gradient of temperatures stretching downwards under the greenhouse and inwards from the base of the walls. Considering the complexity of vertical and horizontal gradients of soil temperature, a simple, generic model is called for. Consequently, we set soil temperature at time t equal to the average outdoors temperature over the last 7 days.

Setpoints

All setpoints are defined to be in effect for a certain time interval, both in terms of calendar period and daily time interval. Thus setpoints may change during the year and during the day.

Proportional setpoints

Some setpoints can take on a range of values depending on an input variable. Here are two examples:

Such proportional setpoints are defined by their range ([100;300] on the y-axis above) in response to an input threshold (15 above) and a threshold band (10 above). The response can be increasing (left figure) or decreasing (right figure).

Proportional setpoints are written in the notation:

S=x[Sx;ΔSx][Sbegin;Send]

The two examples above are, respectively, x[15;10][100;300] and x[15;10][300;100]. This notation might become handy for scientific publication.

Hence, a proportional setpoint involves the definition of four setpoints:

If the band width is set to zero ΔSx=0 then the setpoint can only be in one of two states, either S=Sbegin or S=Send.

Humidity setpoint

The humidity setpoint SHmax (%) and its proportional band ΔSHmax (%) set the upper limit on indoors relative humidity Hin (%). It is used for proportional control of the heating setpoint offset and the ventilation crack setpoint.

This boxscript sets the humidity threshold rhMax to 80 from May to August and otherwise to 90. The proportional band is always 5:

The rhMac/threshold box is a PrioritySignal box, which means that it will check the boxes inside one by one and pick the first one that should be activated according to the current data (maintained by the calendar box; see Boxscript example with calendar). Here as elsewhere, a lonely : denotes left-out lines.

Heating setpoint

The heating setpoint Sheat (°C) is derived from two other setpoints:

A simple boxscript goes to specify this:

Here setpoints/heating/base[value] is 20 from April to September and otherwise 22. The setpoints/heating/huimityOffset[value] is 2 always. The heating setpoint Sheat (°C) will eventually be computed by the heating controller as

Sheat=Sheatbase+Sheathumoffset*

where Sheathumoffset*[0;Sheathumoffset] is the realised humidity offset.

Average heating setpoint

In average climate control, the base heating setpoint Sheatbase (°C) is updated regularly (e.g., at midnight) to reflect the average greenhouse temperature Tin (°C) (e.g., over the last 48 hours). The average heating setpoint Sheatavg (°C) regulates Sheatbase in steps of ΔSheatbase (°C) to keep TinSheatavg:

ifTin>SheatavgthenSheatbaseSheatbaseΔSheatbaseelseSheatbaseSheatbase+ΔSheatbase

This functionality has not yet been implemented.

Ventilation setpoint

The ventilation setpoint Svent (°C) is derived from three other setpoints:

In the boxscript, the humidity offset is given a negative value (here -1.5) to ease the handling of setpoints by the ventilation controller:

The ventilation setpoint Svent (°C) will eventually be computed by the ventilation controller as

(37)Svent=Sheat+SventbaseoffsetSventhumoffset*

where Sventhumoffset*[0;Sventhumoffset] is the realised humidity offset.

Ventilation crack setpoint

At high indoors humidity Hin (%), ventilation must be kept above a certain minimum; this is the 'ventilation crack' setpoint Sventcrack (h1). As its unit implies, it a minimum ventilation rate. However, if the outdoors temperature Tout (°C) is too low then the ventilation crack is not applied.

In this boxscript, the ventilation crack has been set to 0.5. The realised ventilation crack is determined by the ventilation crack controller, governed by the humidity setpoint but also taking into account the ventilationCrack/temperature setpoint. In this example, the ventilation crack will keep its full value (0.5 h1) down to an outdoors temperature of 5+3=2°C. Below 5°C the ventilation crack will be zero.

Screen crack setpoint

A screen crack setpoint would be defined to ensure that screens would not block the free passage of air through the vents, when ventilation were necessary. However, the ventilation model is not detailed enough to correct for the influence of screens on the effectiveness of the vents. Hence, the model does not include a screen crack setpoint.

Growth light setpoints

Growth lights are arranged into banks, which each follow their own logic; growth lights in a bank are switched on/off together. A bank can be operating in any of three modes designated by a code (0, 1 or 10):

The mode is itself a setpoint. Hence, you can change any aspect of the growth light logic through time.

In the following example there are two banks of setpoints, named bank1 and bank2. Banks can be named freely but the same names must be used in growth light controllers and actuators too; furthermore, a bank most hold a child box called mode.

The outcome of these setpoints is for bank1:

For bank2:

Strictly speaking, it is decided by the growth light controller (here the controllers/growthLights/bank1 controller) which variable the thresholds refer to.

Chalk setpoint

The glass may be chalked (whitened) to lessen radiative transmission. The setpoint is a number between 0 and 1 allowing for varying strengths of whitening. The actual strength is set by the chalk controller, which makes it possible to reduce whitening strength with time or due to events such as heavy rainfall. The actual, radiative effect of chalking is determined by the chalk actuator.

The example whitening is fully on from 1 March to 15 April; otherwise it is off:

 

Screen setpoints

Screen setpoints provide setpoint values used by the screen controllers. You should give the setpoint boxes names that will make it easy for you to recall their relations to the screen layers defined in the screen house construction. Please refer to the Boxscript example with cover and screens) which we continue in the boxscript further down.

In general, energy screens will need a threshold below which the curtain should drawn (e.g., in response to outdoors radiation), whereas shade and blackout screens need a threshold above which they will be drawn. Screen setpoints are usually proportional setpoints and thus need a proportional band in addition to the threshold. Both the threshold and the proportional band can be set up with different values for different periods of the year and the day.

Boxscript example with screen setpoints

In this example, the setpoints are divided into three groups, each held inside one box: energy, shade and blackout. A value is set for a certain period by a DateTimeSignal box. For instance, the shade/threshold box holds three DateTimeSignal boxes (lines 20-37). The shade/threshold box itself is a PrioritySignal box, which means that it will check the boxes inside one by one and pick the first one that should be activated according to the current data and time (maintained by the calendar box; see Boxscript example with calendar). The dates beginDate and endDate are inclusive. In this case, the result will be that the output port setpoints/shade/threshold[value] (to spell out the unique path identifying this port) will take on the value 300 in December and January, 400 from 1 February until 15 March and 500 from 16 March until 30 November. Meanwhile, the proportional band at setpoints/shade/band[value] will always be 50.

For easier reading, you can follow my example and always define a PrioritySignal box even though it holds only one DateTimeSignal, which is then a given choice. If none of the boxes held by PrioritySignal matches the current date and time, it will sets its value to zero which is rarely wanted. Hence, it is best always with to have a final DateTimeSignal box that spans the whole year.

The blackout setpoint (lines 46-70) causes the screen to be on from 1 April to 15 June every night from 4 a.m., when outdoor light is above a threshold value of 2. Outside this period of the day and of the year, the threshold value is 9999, i.e. the screen will be off since light (if its in W/m2) will never surpass this value.

CO2 setpoint

The CO2 setpoint sets the desired minimum CO2 concentration (ppm). As an example, the boxscript below setpoint sets CO2 at 1200 ppm in the summer and at 900 ppm for the rest of the year:

At high ventilation rates, CO2 enrichment is a waste. In the boxscript above, the limit has been set at 0.2 h1 with a proportional band of 0.1 h1. Hence between a ventilation rate of 0.2 h1 up to 0.3 h1, the CO2 dispensation rate will dwindle to zero.

The CO2 setpoint is passed on to the CO2 controller, while the ventilation-dependent limitation is applied directly by the CO2 actuator.

Controllers

Controllers add further to the logic imposed by setpoints.

Heating controller

The heating controller takes the base heating set point and adds the humidity offset for heating. Both values are set by the heating setpoint.

In the boxscript below, controllers/heating is a Sum box which adds the base value and the humidity offset (the pipe | operator signifies a union of the matches provide by the two paths setpoints/heating/base[value] and ./humidityOffset[value]).

The offset is calculated by the child box (referred to by a leading period, i.e. ./humidityOffset[value]), which again is a ProportionalSignal governed by the humidity setpoint and the current indoors humidity indoors[rh]. Note how the maximum possible offset is grabbed from the setpoint setpoints/heating/humidityOffset[value]. The offset is increasing with humidity through the threshold band.

Ventilation controller

The ventilation controller works very much like the heating controller, except it adds up three values (cf. 37):

The three values are taken from three different sources: the heating controller, the ventilation setpoint and the realised humidity offset computed by its own child box. The ventilation actuator must take both the ventilation controller and the ventilation crack controller into account when determining the actual ventilation rate.

Ventilation crack controller

The ventilationCrack controller first restricts the size of the crack by the temperature limitation imposed by outdoors temperature (lines 9-15) (in BoxScript, children are updated before parents). This crack value ./temperatureLimited[value] is then moderated further by way of the humidity setpoint to arrive at the final crack controllers/ventilationCrack[value].

 

Growth light controllers

There may be several growth light controllers, one for each bank of growth lights. A growth light controller switches the growth lights in the bank on/off according to the growth light setpoints. To continue the setpoint example, we can define the following two controllers:

Due to the complicated control logic, a dedicated GrowthLightController is needed for each bank. The bank1 controller refers to the bank1 setpoints (lines 5, 7 and 8) and likewise for bank2 (line 12). Only bank1 needs the inputs thresholdValue, thresholdLow and thresholdHigh; bank2 is not threshold-controlled. The optional minPeriodOn can be set to obtain a certain minimum periods (in minutes) for the lights being switched on.

CO2 controller

The CO2 controller is a PID controller which controls the CO2 injection system to regulate the CO2 towards the CO2 setpoint.

A PidController box has three parameters:

These are used to steer the sensedValue towards the desiredValue. This boxscript provides an example in which only proportional control Kprop is applied:

PidController produces the controlVariable output which tends towards zero, when the P+I+D terms tend towards a zero sum. Here, it used as an indicator how much the CO2 injection rate should change (line 4). The co2[change] output is used as an input by the CO2 actuator.

Screen controllers

The screen controllers use the screen setpoints together with inputs provided, e.g., from the outdoors box to compute desired screen states between 0 (fully off) to 1 (fully on). As mentioned earlier, there is no screen crack setpoint to take into account.

Boxscript example with screen controllers

With the definitions given in the boxscript example with screen setpoints, we could write the following three controllers, which all responds to outdoors radiation (W/m2) provided by outdoors[radiation]:

Note how the increasingSignal input tells whether the signal computed by a ProportionalSignal box is increasing or decreasing in response to the input. There is also a minSignal input but this defaults to 0, which is just what we want.

All boxes of the Signal classes (ProportionalSignal, PrioritySignal, etc.) provide their about under two equivalent names, value and signal.

Chalk controller

The chalk controller may provide mechanisms to modulate the whitening strength set by the chalk setpoint. Most simply it just passes on the setpoint:

The chalk actuator retrieves the current whitening strength from /controllers/chalk[value].

Actuators

Growth light actuators

Each bank of growth lights is controlled by a growth light controller for which there exists a growth light actuator with the same name (e.g., bank1, bank2)

A growth light actuator is described by the parameters:

The efficiency of lamps will degrade with their age. The product ϵlampoptϵlampage expresses the actual efficiency.

See the Growth light layer section for details on growth light function. Then outputs of all sets of growth lights are added and considered one layer.

 

Screen actuators

Each screen actuator is governed by a controller. It has a lagPeriod input, which is the time taken (in minutes) to drag the screen all the way from position 0 to 1 (or from 1 to 0). If the lag period is much smaller than the simulation time step, it is better to set it to zero to avoid unrealistic delays in effectuating screen settings.

Boxscript example with screen actuators

To continue the boxscript example with screen controllers, we could write

 

Heat pipe actuators

There can be several sets of of heat pipes installed. Each set of heat pipes is described by the parameters:

See the Heat pipe layer section for details on heat pipe function. Then outputs of all sets of heat pipes are added and considered one layer.

The actual pipe inlet temperature Tpipeinlet (°C) is governed by the energy budget model according to the relevant controllers.

Ventilation actuator

Vents are characterised by their width and length and the number of vents installed. If vents of different dimensions are installed, they are specified as different sets. The total area of vents Avent (m2) is used to calculate the maximum ventilation rate vmax (h1) (eq. 17). Due to the simplicity of the ventilation model, no specification is needed for the position of the vents.

The actual ventilation rate v (h1) (see advection) is governed by the energy budget model according to the relevant controllers.

CO2 actuator

The CO2 actuator injects CO2 at a rate of cCO2 (g/m2/h). The actuator is limited by its maximum injection capacity CCO2 (g/m2/h). In the boxscript below it has been set to 60:

When deciphering this boxscript, remember that boxscripts are interpreted child-first. Hence injectionRate in lines 9-15 comes first. Here, the maximum injection rate is adjusted for the current ventilation rate, which might limit it below the set maximum of 60 (line 13). This value is taken by the parent box (lines 5-16) as an input (line 8) to set the maximum ventilation rate arrived at by PID control. Accumulator is a box that in every time step adds the change input (line 6) to its current value. The value can, however, never exceed the bounds set by minValue and maxValue (lines 7-8). The change input was computed by the CO2 controller by PID control. The resulting valueis finally set as the injectionRate in line 4. Thus the CO2 budget model can retrieve the current CO2 injection rate from the port actuators/co2[injectionRate].

Chalk actuator

The current strength of whitening (ϵchalk[0;1]) is provided by the chalk controller. In addition, the chalk actuator takes two parameters for the additional reflectivity provided by the whitening, ρchalksw[0;1] and ρchalklw[0;1] for short-wave and long-wave radiation, respectively. The corrections, ϵchalkswρchalk and ϵchalkρchalklw, are applied to the cover of all six greenhouse faces for both incoming and outgoing radiation.

As an example, if the cover has the properties ρsw=0.1, τsw=0.8 and αsw=0.1, and the chalk effect is ϵchalkρchalksw=0.2 then the corrected properties of the cover for incoming radiation will be ρsw=0.1+0.2=0.3, τsw=0.80.2=0.6 and αsw=0.1. Likewise, for the parameters for outgoing radiation. The actuator enforces limits to ensure that the corrected properties are all between 0 and 1, and that they sum up to 1.

Here is a boxscript example:

Budgets

Energy budget

We build the energy budget by adding the processes of heat exchange one by one. The step-wise construction is reflected in R scripts that progressively includes more and more of the energy budget.

Sources of radiation

Radiation enters the model as layers may emit radiation downwards from the bottom (Eisw,Eipar,Eilw) and upwards from the top (Eisw,Eipar,Eilw). While only some layers act as a primary source of short-wave radiation radiation (sky and growth lights), all layers may transmit and reflect short-wave radiation received from neighbouring layers. In contrast, all layers are primary sources of long-wave radiation. For most layers, the long-waved radiation is determined by the Stefan-Boltzmann Law (eq. 12) ('S-B Law' in the table below).$

The primary emission plus the transmitted and reflected radiation defines the net flows of radiation downwards (Fisw,Fipar,Filw) and upwards (Fisw,Fipar,Filw) from each layer. For the sky, all upwards radiation is zero Eir=0; likewise for the floor, all downwards radiation is zero Eir=0.

LayerTemperature
(Ti)
Short-waved radiation
(Eisw,Eisw,Eipar,Eipar)
Long-waved radiation
(Eilw,Eilw)
Skydriving variabledriving variablef(Ti) (S-B Law)
Coverstate variable, f(H˙i)0f(Ti) (S-B Law)
Screenstate variable, f(H˙i)0f(Ti) (S-B Law)
Growth lightn.a.state variable, f()state variable, f()
Plant canopystate variable, f()0f(Ti) (S-B Law)
Heat pipestate variable, f()0state variable, f()
Floorstate variable, f(H˙i)0f(Ti) (S-B Law)

For the sky both its temperature and short-wave emission are driving variables (usually read from a weather log file). Its downward long-wave emission is calculated from its temperature.

For cover, screen and floor layers, temperature is updated according to the heat balance H˙i (W/m2) and the heat capacity Ci (J/K/m2),

ΔTi=H˙iΔtCi

where Δt (s) is the simulation time step and ΔTi (K) is the change in layer temperature.

For growth light, temperature is not tracked. A dedicated sub-model, denoted as a function with numerous inputs f(), computes the emission of radiation (short-waved, PAR and long-waved). Dedicated sub-models also calculated the temperature of the canopy and heat pipe layers, as well as the long-waved radiation emitted from heat pipes.

Short-wave radiation

As radiation is emitted, absorbed, reflected and transmitted by the layers, a net radiation flux will result downwards and upwards for each layer i. The following algorithm, here applied to short-wave radiation, will also be used for long-wave radiation and PAR.

The net radiation fluxes throughout the stack of layers are resolved by resolving them two layers at a time. Consider layer i on top of the layer i+1 below:

energy-budget-layers

Our first aim is to find out, how much of the downwards radiation from the upper layer Fi (W/m2), e.g Fisw, is absorbed by the lower layer? If we denote the net absorption from above α^i+1, we immediately have

α^i+1=αi+1

However, a fraction of the light will be reflected by the upperside of the lower layer (ρi+1) and then re-reflected from the underside of the upper layer (ρi). We must add that bit,

α^i+1=αi+1+αi+1ρi+1ρi

and so forth as photons keep bouncing between the two layers,

α^i+1=αi+1+αi+1(ρi+1ρi)1+αi+1(ρi+1ρi)2+=αi+1j=0(ρi+1ρi)j

Since the reflectivies are less than one, the infinite series will converge as

ri=j=0(ρi+1ρi)j=11ρi+1ρi

and we get the amount absorbed by the lower layer Ai+1 (W/m2),

Ai+1=α^i+1Fi=riαi+1Fi

The net transmission through the lower layer τ^i+1 follows the same logic and we arrive at

τ^i+1=riτi+1

On every bounce of the radiation, fractions will be lost to absorption by and transmission through the upper layer. Thus the calculations for the underside of the upper layer are the same, except that the we start out with what's first reflected upwards from the lower layer (ρi+1). Hence, we get

α^i=ρi+1riαiτ^i=ρi+1riτi

One consequence of these reflections is that part of the downwards radiation from the upper layer Fi will be absorbed by the layer itself from below Ai (W/m2):

Ai=α^iFi=ρi+1riαi+1Fi

Let's check that we have now accounted for the fate all radiation downwards from layer i. We have split it into radiation absorbed by and transmitted through either of the to layers. In total we get

α^i+1+τ^i+1+α^i+τ^i=αi+1+τi+1+ρi+1(αi+τi)1ρi+1ρi=(1ρi+1)+ρi+1(1ρi)1ρi+1ρi=1ρi+1+ρi+1ρi+1ρi1ρi+1ρi=1

Thus we have accounted for it all.

Having accounted for the fate of the downwards radiation from layer i, we can now repeat the calculation for the next pair of layers, i+1 and i+2, and so forth until the bottom of the stack. The downwards flow from any layer i is the inherent emission Ei (i.e. for a layer of growth lights) plus the transmission from the layer above τ^iFi1,

Fi=Ei+τ^iFi1

After we have finished the distribution of downwards radiation, we repeat the procedure this time going from the bottom to the top following the upwards radiation,

Fi=Ei+τ^iFi+1

We repeat the calculations downwards-upwards, on each pass accumulating in each layer the radiation absorbed from above Ai and below Ai. The calculation stops when nearly all radiation has been absorbed, e.g. until iFi+iFi<106 W/m2. Note that the inherent emission downwards Ei and upwards Ei is only added to the respective flows, Fi and Fi, on the first pass through the calculations.

PAR

The PAR budget is resolved just as the short-wave budget above. However, PAR is not part of the energy budget. Its only role is to provide energy for photosynthesis in the canopy layer. It needs a special treatment because the emission of PAR (Eipar,Eipar) and short-wave radiation (Eisw,Eisw) will differ between layers, such us sky and growth lights.

Computation example

The 1-energy-budget-sw.R script demonstrates the calculations. Here the light layer emits 100 W/m2 downwards (or 100 μmol/m2/s; the calculations are the exact same). The zero emission of short-wave radiation from the sky indicates a night situation. Greek letters are represented by their Latin equivalent and primes have been replaced by underscores:

The radiative parameters for the plant layer was computed for Lai=1.9 and k=0.65 giving α=0.71, ρ=0.05 and τ=0.24 (eqs. 24 to 26).

After the first distribution of radiation downwards, we have

We can see that 71.0+14.4=85.4% of the radiation has been absorbed by plants and floor together, while 5.0+9.6=14.6% has been reflected upwards by the same layers. We carry on distributing these 14.6% upwards and get

This led to a total of 2.0+0.1+0.7+6.8=9.6% (A_ column) being absorbed on the underside of various layers (sky, glass, screen, plant). Meanwhile, reflection downwards led to the remaining radiation of 0.1+4.4+0.5=5.0% (F column) hitting som layers (screen, light, heating) from above.

If we set the precision, so that the calculations down and up will be repeated until Fi+Fi<106 W/m2, it turns out that seven iterations are needed, and we get the final distribution,

Thus 74.3+7.3=81.6% of the radiation ended up in the plant canopy (which is the intended target of growth lights). Most of the remaining light was absorbed by the floor (15.4%), while only little ended up in glass and screen. Only 2.2% escaped to the sky; this is the light that you can see passing a greenhouse at night.

Long-wave radiation

The same algorithm as for short-wave radiation above is used to resolve the fate of long-wave radiation, i.e. its distributed among the layers. However, the absorption of radiation (both short-wave and long-wave radiation) will cause a change in the temperature of the layer Ti (°C) according to its heat capacity Ci (J/K/m2). This is taken into account by an additional step in the algoritm, which updates layer temperature as radiation is absorbed.

Note that the distribution of short-wave radiation remains constant because it is temperature-independent. Only the long-wave emission of the layers will change and converge, as a solution is found in which all temperatures are stable and by that, the heat exchange by radiation is in a steady state.

Computation example

The 2-energy-budget-lw.R script demonstrates the calculations, which follow the same logic as for short-wave radiation above but with some extra parameters and calculation steps.

Step 1. Initial set-up

Here is the initial set-up. Columns T and C have been added to hold temperature and heat capacity. Incidentally, Cglass was set to 8400 J/m2/K by a mistake. The intention was to be consistent with the example in eq. 18, where it is found to be 10416 J/m2/K.

We have transferred the absorbed short-wave radiation from the example above (columns A and A_). The emissions (E and E_) have been set according to layer temperatures and emissivities, applying Stefan-Boltzmann's Law (eq. 12), except for light and heating layers. The long-wave emission from light (15 W/m2 both up and down) and heating (75 W/m2 both up and down) would in the full model be calculated by specific sub-models (see growth light and heating). Note that we use absorptivity as a stand-in for emissivity (see layers).

The radiative parameters for the plant layer was computed for Lai=1.9 and k=1.0 (leaves are impenetrable ('black') to long-wave radiation) giving α=0.79, ρ=0.01 and τ=0.20 (eqs. 24 to 26). Glass, likewise transmit no long-wave radiation τ=0.

It worth noting that screens may have an asymmetric emissivity for long-wave radiation (ϵilwϵilw). Here we used the radiative parameters exemplified earlier (see screen layer); however, Cscreen=2280 J/m2/K indicates a sturdier material than the example shown in eq. 23. A sky emissivity of 1 (αsky=1) suggests a sky with no cloud cover.

The temperature of layers that have been given an infinite heat capacity, does not change in temperature by the absorbed radiation. However, the plant temperature is kept fixed in this budget for now, until we have included the full plant canopy sub-model (see canopy temperature), which among other things calculates leaf temperature.

Step 2. Absorb radiation

We use the same algorithm as for short-wave radiation to distribute the flows of long-wave radiation. We arrive at this solution of the long-wave plus short-wave radiation budget:

We left step 1 with Ai+Ai=100 W/m2 (columns A and A_) and Ei+Ei=2324.2 W/m2 (columns E and E_). After step 2 the long-wave emission has been added, and we have Ai+Ai=2424.2 W/m2.

Step 3. Update temperature

The net absorbed radiation for any layer Airad is found by subtracting the emission from the absorption. For example, we have for the glass layer (i=2),

A2rad=A2+A2E2E2=239.3+345.2344.0344.0=103.5 W/m2

How much of a change in glass temperature this will cause, depends on the time step. Let's say the simulation time step for the whole model is Δt=180 s but that we sub-divide this into n smaller time steps, when computing the energy budget to improve the numerical precision. With n=6 we achieve a change in glass temperature ΔT2 (K) of

ΔT2=A2radC2Δtn=103.5 W/m28400 J/K/m2180 s6=0.37 K

In general, we should set n small enough to achieve an acceptable, small |ΔTi| for all i. To obtain a fast execution time, we can adjust n in every simulation time step as needed. If we want a maximum temperature change of, say, 0.5 K in all layers, we can find n by

n=max(|ΔTi|)0.5 K

where ΔTi are the temperature changes calculated with n=1. The ceiling operators ... round up to the nearest integer. This is how n=6 was determined for this example.

We update the temperature of screen and floor in the same manner and get

In summary, the glass was cooled by 0.37 K, while screen and floor was heated by 0.48 K and 0.08 K, respectively, during the 30 s time step. The net absorption by the plant canopy is ignored for now. The absorption columns A and A_ have been set to zero too reflect that the absorbed energy has been transformed to heat.

Step 4. Update long-wave emission

Since the temperature of glass, screen and floor has changed, we must update their long_wave emission in columns E and E_:

Step 5. Keep invariant emissions

The short-wave absorption is unaltered and is once again (as in step 1) entered into columns A and A_. Likewise, the long-wave radiation of light and heating remain unchanged in columns E and E_:

Following steps

After step 5 we return to step 2 and loop through the calculations, step 2 to 5, for a total of n=6 iterations to achieve the final result after Δt=180 s:

Verification

The total net absorption for any layer i over the the j=1..n minor time steps ΔAinet (J/m2) is

(38)ΔAinet=Δtnj=1nAijsw+Aijsw+Aijlw+AijlwEijlwEijlw

where we for convenience have left out the emission of short-wave radiation, which amounts to

Eijsw+Eijsw=100 W/m2

The energy accumulated by the sky and plant layers were A1net=18471 J/m2 and A5net=25690 J/m2 (not shown in the table above), respectively. To verify that the algorithm is correct, we can calculate the net absorption for the glass, screen and floor layers based on their temperature change ΔTi during the whole simulation step. We get

ΔA2net=ΔT2C2=(12.915.0) K8400 J/K/m2=17726 J/m2ΔA3net=ΔT3C3=(20.518.0) K2280 J/K/m2=5633 J/m2ΔA7net=ΔT7C7=(12.512.0) K42000 J/K/m2=18332 J/m2

Noting that the light and heating layers absorb no radiation, we get

ΔA4net=(E4+E4)Δt=30 W/m2180 s=5400 J/m2ΔA6net=(E6+E6)Δt=150 W/m2180 s=27000 J/m2

We can now calculate the energy balance of the whole system,

iΔAinet=1847117726+56335400+2569027000+18332=18000 J/m2

This can be compared against the total energy put into the system, namely the 100 W/m2 short-wave radiation. Over the the whole simulation time step this gives

100 W/m2180 s=18000 J/m2

which balances the whole budget to zero and thereby verifies the precision of our calculations.

Convective and conductive heat transfer

Unlike the exchange of radiation, which occurs between layers, convective heat transfer occurs between a layer and its surrounding volume of air. The cover exchanges heat with the outdoors air on the outer surface and the inside air on the inner surface. The floor exchanges convective heat with the indoors air. All other layers exchange convective heat with the indoors air on both upper and lower surfaces. The floor in addition exchanges conductive heat with the outdoors soil.

All convective/conductive heat transfers are defined by the U-value of the layer (eq. 13). For the outer surface the U-value depends on the wind (eq. 22), while for all the inner surfaces, U is set to 1.2 W/K/m2. We assume a good insulation of the floor against the soil setting Ufloor=0.1 W/K/m2.

We will keep track of the convective/conductive heat absorbed by each volume (Aoutconv,Ainconv,Asoilconv) (W/m2) to check that the total energy budget is consistent, i.e. the energy of the whole system is conserved. To update the indoors temperature Tin by the absorbed heat Ainconv, we apply the area-specific heat capacity of the air,

Cair=1020JKkg1.19kgm33.94m3m2=4780JKm2

where 3.94 m is the average indoors height (eq. 10). The change in indoors temperature ΔTin (K) over a time step Δt (s) is

ΔTin=AinconvCairΔt

Computation example

The convective/conductive heat transfer between layers and volumnes is added to the short-wave and long-wave radiation budget, as demonstrated in the 3-energy-budget-convection.R script.

Step 1. Initial set-up

We have added Ui and Ui (columns U and U_) and heat absorbed by convection/conduction Aiconv and v (columns H and H_) to the layer budget. In addition, we need a separate budget for the three volumes with columns T, C and H representing Tj, Cj and Ajconv, respectively, where j subscript outdoors, indoors and soil:

All the functionality of the short-wave (energy-budget-sw.R) and long-wave budgets (energy-budget-lw.R) is kept. The heat convective/conductive heat exchange between layers and volumes can simply be been added.

The initial set-up now includes the convective heat emitted by growth lights and heat pipes. This has been set to Hlight=Hlight=5 W/m2 for growth lights and Hheat=Hheat=50 W/m2 for heat pipes, both downwards and upwards. This leads to a total gain of the indoors air of Tin=110 W/m2. The plant layer is given a provisional U-value of 0.1 W/K/m2 (that is, until we have added the plant canopy sub-model to the energy budget).

Step 2. Absorb radiation

This step only involves the layers and gives the same result as step 2 in the long-wave radiation budget; the E and E_ columns are distributed into the A and A_ columns:

Step 3. Update heat fluxes

The convective/conductive heat fluxes are governed by the layer U-values. This step affects both layers and volumes, except the light and heating layers, for which the heat fluxes are fixed:

We can check that all heat fluxes (columns H and H_) add up to zero:

Layers:i=17Aiconv+Aiconv=79.8 W/m2Volumes:Aoutconv+Ainconv+Asoilconv=79.8 W/m2

In other words, the layers are giving off heat to all three volumes, most of it to the indoors and the outdoors air. Only little disappears through the well-insulated floor.

Step 4. Update temperature by radiation

This step corresponds to step 3 for in the the long-wave radiation budget. We get the same resulting layer temperatures as before, as the A and A_ columns are converted to temperature changes and added to the T column:

Step 5. Update temperature by heat fluxes

This step is similar to step 4, but now it is the H and H_ columns that are converted to temperature changes and added to the T column. This step involves both the layers and the indoors volume. Of the layers, only the screen layer changed discernably:

Step 6. Update long-wave emission

Since the temperature of glass, screen and floor has changed, we must update their long_wave emission in columns E and E_:

 

Step 7. Keep invariant fluxes

The short-wave absorption is unaltered and is entered into columns A and A_, likewise, the heating fluxes in columns H and H_:

 

Following steps

After step 7 we return to step 2 and loop through the calculations, step 2 to 7, for a total of n=6 iterations to achieve the final result after Δt=180 s:

The columns SumNetAH and sumNetH corresponds to the net absorption ΔAinet defined for layers and volumes, respectively (eqs. 39 and 40).

When considering the loss to the outdoors, it should be noted that ΔAskynet=18484 J/m2 is the radiant heat loss, while ΔAoutnet=4741 J/m2 is the heat lost by convection.

Verification

The total net absorption for any layer i over the the j=1..n minor time steps ΔAinet (J/m2) is an extended version of eq. 38,

(39)Layers:ΔAinet=Δtnj=1nAijsw+Aijsw+Aijlw+AijlwEijlwEijlw+Aijconv+Aijconv

where again we for convenience have left out the emission of short-wave radiation, Eijsw+Eijsw=100 W/m2.

Similarly for the three volumes (i=1..3), which only exchanges heat through convection/conduction:

(40)Volumes:ΔAinet=Δtnj=1nAijcon

We get a total balance of

Ainet[layers]+Ainet[volumes]=5483+12517=18000 J/m2

which again matches the emitted short-wave energy, 100 W/m2180 s=18000 J/m2.

Advective heat transfer

Advective heat transfer is effectuated by ventilation. This causes a transfer of energy, in cold climates usually a loss, from the indoors air to the outdoors. This is due to a loss of both sensible heat (the air being colder outside than the inside) and of latent heat (the air being drier outside than inside). The advective heat transfer treated here concerns only the sensible heat while the transfer of latent heat is treated further down (after the water budget has been worked out).

Computation example

The 4-energy-budget-advection.R script continues the example from above with a simulation time step Δt=180 s, broken into n=6 minor time steps and with temperature indoors at an initial Tin=25°C and outdoors fixed at Tin=10°C. We choose a ventilation rate of 3 h1. Over a minor time step of 30 s, this leads to a cooling of the greenhouse air ΔTin=0.37°C (eq. 15) corresponding to an advective heat flux of Ainadv=59.0 W/m2 (eq. 16).

Our initial set-up for the budget has been extended with a V column for the volumes to hold Aiadv:

The steps to resolve the energy budget follows those above with the addition of the ventilation flux (column V) in every minor time step. Thus we arrive at this solution for the budget:

Beforehand, the budget resulted in an increase in Tin due to a net transfer of heat from the layers to the greenhouse air. Now, with ventilation included Tin is lowered due to the increased heat loss to the outdoors air.

The total net absorption of layers is still described by eq. 39, while for volumes we must add (to eq. 40) the contribution of ventilation:

(41)Volumes:ΔAinet=Δtnj=1nAijconv+Aijadv

Again we achieve the expected energy balance,

Ainet[layers]+Ainet[volumes]=4162+13838=18000 J/m2

Canopy temperature

So far, we have set the U-value of the canopy to 0.1 W/J/m2 and fixed its temperature at 24°C. A more correct model of canopy sets its U-value to zero and updates its temperature by eq. 27, which we will do now.

Computation example

The code needed only little change to implement these changes. This is the outcome of the resulting 5-energy-budget-plant.R script:

The relative humidity was fixed at 90% for these calculations. It can be seen that the plant temperature, previously fixed at 24°C, has risen to 25.4°C. This again has caused the screen temperature to rise from the previous 21.4° to now 21.6°. The temperature of other layers changed less.

The total energy budget changed only little and is still in balance:

Ainet[layers]+Ainet[volumes]=3793+14207=18000 J/m2

Water budget

The dominating source of water vapour in the greenhouse is canopy transpiration. The main sinks of water vapour are condensation and ventilation. Condensation on the glass is welcome as it lowers air humidity and thereby lowers the risk of condensation on the plant as well. Condensation on the plants is unwanted as it may lead to disease outbreaks. We consider water condensing on the plants an exception and include only condensation on the glass in the model.

The water film forming on the glass is very thin, only about 16 μm (Pieters el al., 1996). During prolonged condensation (typically during the night) water will run off the glass and end up in the gutter. Evaporation of the water film thus makes only a tiny contribution to the water budget and is left out of the model.

The water budget model consists of three components: transpiration from the canopy, condensation on the inside of the glass, and mass transport of water between indoors and outdoors air.

The figures in this sub-section were generated by the 5-water-budget.R script.

Transpiration

TranspirationHplanttrans (kg/m2/s) is modelled by eq. 33. It is highly dependent on the radiation absorbed by the canopy, as seen in the figure,

Condensation

Condensation occurs when the glass temperature falls below the dew point of the indoors air. The saturated water vapour pressure at the temperature of the glass is Hsatpres(Tglass) according to eq. 28, while the water vapour pressure of the indoors air is

Hinpres=Hinrel100%Hsat(Tin)

The absolute humidity Habs (kg/m3) can be calculated from pressure and temperature using the Ideal Gas Law:

(42)nV=PRT

In our terms,

(43)Habs=MwHpresR(T+T0)

where Mw=18 g/mol is the molar mass of water, R=8.314 m3 Pa/mol/K and T0=273.15 K. The relation looks like this

It is seen how warmer air can hold more water vapour.

The saturated absolute humidity only depends on temperature (eqs. 28 and 43),

Hsatabs(T)=MwHsat(T)R(T+T0)

The condensation rate at the glass Cglasscond (kg/m2 ground/s) is proportional to the difference between the saturated absolute humidity at the glass and the absolute humidity of the air (GCC, eqs. 3.4.18, 3.4.23),

(44)Cglasscond=CaicHinabsHsatabs(Tglass)0

with units

[kgm2 grounds]=[m2 coverm2 ground][m3m2 covers][kgm3]

where Cai is the cover/ground area index (eq. 9) and c (m/s) is an empirical parameter. For condensation on screens, Cai must be replaced with the ratio of the screens-drawn area to ground area.

Pieters et al. (1996) estimated c from data collected in a controlled experiment with Tin=20°C and Tout=10°C at relative humidity Hinrel=(70%,80%,85%,90%). The temperature of the glass was not measured. We transformed the experimental Hinrel to Hinabs and tentatively assumed a glass temperature between the inner and outer temperature, Tglass=(13,15,17)°C. We estimated c=2103 m/s, as seen in the figure, where the points show the experimental data:

This estimate of c is in accord with that in GCC (eq 3.4.23), which has Caic=3103 m/s. With a typical value of Cai=1.5 m2/m2, the two estimates are equivalent.

This figure shows the condensation rate, depending on glass temperature, at indoors temperature Tin=24°C and two different indoors humidities Hinrel=(60%,90%),

The vertical line at 24°C highlights that at high indoors humidity (90%), just a slightly cooler glass will accumulate condensation. The condensation rate is in the same range as transpiration, e.g., at 90% humidity and a glass temperature of 10°C, the condendation rate is 0.02 g/m2/s (read off the figure) giving 1.2 g/m2/min.

 

Ventilation

Ventilation causes mass movement of water vapour between outdoors and indoors air. To calculate the resultant change in indoors absolute humidity ΔHinabs (kg/m3) during one time step Δt (s), eq. 15 can be rewritten as

(45)ΔHinabs=[HoutabsHinabs(0)][1exp(v3600 s/hΔt)]

where Hinabs(0) (kg/m3) is the indoors absolute humidity at the beginning of the time step.

The average rate of water lost by mass movement during the time step is

(46)ΔHinabsVaiΔt[kgm2s]

where Vai is the volumetric area index (eq. 10).

This is depicted in the following figure at a ventilation rate of v=3 h1 and with temperatures inside Tin=24°C and outside Tout=10°C. The relation is shown at two different indoors humidities and varing outdoors humidity,

The cooler outdoors air can hold only little water. Hence, ventilation is an effective means of lowering indoors humidity. The effect is larger with more humid indoors. The flux of water vapour is negative because water is lost from the greenhouse. Since the humid indoors air holds latent heat originally captured by transpiration, ventilation also functions as an indirect route of venting excess heat.

Latent heat

The latent heat carried by water vapour is released during condensation. Knowing the condensation rate at the glass Cglasscond (kg/m2 ground/s) (eq. 44), it is straightforward to calculate the heat, which is assumed to be absorbed by the glass immediately at a rate of Aglasscond (W/m2 ground),

Aglasscond=λCglasscond

The heat of vapourisation λ=2454 kJ/kg determines the energy released by the phase shift.

Computation example

The model was extended with a water budget shown as three columns (in g/m2) in the output, Tr (transpration), Cn (condensation) and Mv (mass movement). This is the outcome of the final 7-energy-budget-complete.R script:

 

The output shows the outcome after one time step Δt=180 s, still divided into n=6 minor time steps for the calculations. Thus, 8.0 g/m2 (Tr column) was transpired by the crop. Meanwhile ventilation got rid of 2.5 g/m2 (Mv column) of the the transpired water.

The condensed water amounted to 0.00440 g/m2 (Cn column). Multiplication by the latent heat of evaporation λ=2454 J/g leads to the amount of released energy,

0.00440 g/m22454 J/g=10.8 J/m2

which appears in the SumCn column for the glass layer. The increase in glass temperature caused by condensation alone, taking the heat capacity of the glass layers Cglass=8400 J/K/m2 into consideration, amounted to very little,

ΔTglass[condensation]=10.8 J/m28400 J/K/m2=0.0013 K

Indoors humidity was initialised high at 90% but increased even more to 98.1%, despite the lower outdoors humidity at 60%. In a real greenhouse, ventilation would need to be increased to push humidity down.

The total energy budget changed very little and remained in balance:

Ainet[layers]+Ainet[volumes]=3793+14207=18000 J/m2

It is worth noting that the latent heat, which escaped by ventilation, is not part of the energy budget. In this example it amounted to

2.532 g/m22454 J/g=6214 J/m2

These 6214 J/m2 were originally captured as radiative heat by the canopy and turned into latent heat. They do not constitute an extra loss of energy from the greenhouse, i.e.they do not act as an extra cooling mechanism, as humid air leaves through the vents.

CO2 budget

We use the Ideal Gas Law earlier on (eq. 42) to calculate condensation. Here we apply it to convert CO2 concentration from [CO2]ppm (ppm) to absolute concentration [CO2]abs (g/m3):

[CO2]abs=[CO2]ppm106PRTMCO2=[CO2]ppm106101325 Pa8.314 m3 Pa/K/mol44 g/mol=1.829mg/m3ppm CO2[CO2]ppm

As an example, with greenhouse dimensions of 10000 m2 area and a volume of 39384 m3, how much CO2 would we need to increase indoors CO2 from 400 ppm to 900 ppm? This would amount to

1.829mg/m3ppm CO239384 m3500 ppm CO2=36.0 kg CO2

or

36.0 kg CO210000 m2=3.60 g CO2/m2

or

36.0 kg CO239384 m3=0.915 g CO2/m3

However, the plants will take up CO2 which must be compensated for to keep up the CO2 concentration. For example, a tomato crop grown at a leaf area index of 3.0 with a CO2 uptake rate of 1.29 mg CO2/m2 leaf/s (Heuvelink, 1995) will need

1.29 mg CO2m2 leafs3m2 leafm2 floor3600shg1000 mg=13.9 g CO2/m2/h

or

13.9 g CO2/m2/h10000 m239384 m3=3.54 g CO2/m3/h

Meanwhile, CO2 is lost to the outdoors even with all windows closed due to leakage, maybe at a rate of some 0.5 to 1.0 per hour, depending on wind conditions.

The whole CO2 budget can be expressed in the differential equation,

d[CO2]inabsdt=c+v([CO2]outabs[CO2]inabs)

with

Integration gives us the following solution of the [CO2]inabs time trend with fixed rates c and v,

(47)[CO2]inabs(t)=exp(vt){[CO2]inabs(0)+(cv+[CO2]outabs)(exp(vt)1)},v>0

The figure below illustrates eq. 47 with the following parameter settings:

The time course of the indoors CO2 is shown on the y-axis as [CO2]inppm, i.e. converted to ppm:

Note that c=0 implies that CO2 injection and fixation cancel out. In other words, the injection rate equals the fixation rate (3.54 g CO2/m3/h in the example above). Hence c=0.5 implies an injection rate of 4.04 g CO2/m3/h. With no net gain in CO2 (red curves) the indoors concentration falls from the initial 500 ppm to the outdoors 400 ppm, whereas with a net gain (blue curves) the indoors concentration approaches an asymptote. With increased ventilation (dashed curves) the blue asymptote is lowered, whereas the red asymptote is approached more quickly.

In the full model the CO2 setpoint, controller and actuator will collaborate to adjust the injection rate, so that the indoors CO2 concentration is maintained at the setpoint. The injection rate cinj (g/m2/h) (see CO2 actuator) and the fixation rate Pncanopy (μmol /m2 ground/s) (see Canopy photosynthesis). The exchange rates between units are

1[gm2h]=1[gm2h]AghVgh[m2m3]=AghVgh[gm3h]

and

1[μmolm2s]=1[μmolm2s]AghVgh[m2m3]44.01[gmol]106[molμmol]3600[sh]=0.1584AghVgh[gm3h]

with greenhouse area Agh and volume Vgh.

Implementation notes

Optimisation

The combined energy and water budget is a central part of the model. It is where most computation takes place and where most execution time is spent. Therefore it is worth optimising its performance. This can be achieved by varying the number of minor time steps n, which determines the effective time step Δt/n of the budget computation. One possible algorithm for determining n is

  1. Set n=1.

  2. Compute the budget for one minor time step of length Δt/n and determine the maximum change in temperature of any layer or volume ΔTmax.

  3. If the max. temperature change exceeds a critical limit ΔTmax>ΔTcrit then re-calculate n by

    n=ΔTmaxΔTaccept

    where ΔTaccept is the acceptable change and .. rounds up, Then go back to step 2.

  4. Complete the budget computation for the remaining 2..n minor time steps and determine ΔTmax for the whole procedure.

  5. Re-calculate n by the same equation as in step 3 and go back to step 2. This completes one simulation time step of length Δt for the budget.

This algorith will adjust n both upwards and downwards as deemed sufficient by the limits, which could be for example ΔTaccept=0.5 K and ΔTcrit=1.5 K. The backtracking that could happen in step 3 should be avoided as much as possible to save time. The expectation is that during the night n would be 1. In the daytime it would be higher and when abrupt changes occur in the driving variables (e.g., when growth lights are switched on), it would be even larger for a short while. The simulation time step Δt could be rather large, maybe 5 or even 10 minutes. The proper tradeoff between precision and execution time must be found by simulation.

Photosynthesis

Leaf photosynthesis

The photosynthesis model is identical to the implementation in the Photosyn function found in the plantecophys R package (Duursma, 2015). The Photosyn function was re-coded in C++ keeping only the functionality needed for our purpose. However, boundary layer conductivity gb (mol/m2/s) was added to the function, which originally accounted only for stomatal conduction gs (mol/m2/s). Assuming that boundary layer and stomata work like resistors in a series, they are additive,

1gt=1gb+1gsgt=gbgigb+gi

where gt (mol/m2/s) is the total conductance. All conductances are for water vapour; they are translated into CO2 conductances where needed in the model.

In the overview of variables and parameters below, the mathematical symbols used in this treatise (e.g., Iabs) are listed together with their name in the code (e.g., PPFD). The re-coded versions of the original Photosyn functions can be found in the photosynthesis.R script and, equivalently, in the photosynthesis.cpp source code. All m2 units are per m2 leaf. This is stated explicitly only where there is room for confusion.

The photosynthesis.R script produces this figure to illustrate the output of the photosynthesis function (parameter values can be found in the code):

We see the expected form of the light response curve with positive effects of light, temperature and CO2. Internally, the model lets the more limiting process, electron transport or carboxylation rate (as judged from Pj and Pc), decide the net outcome Pn.

Canopy photosynthesis

The energy budget model computes the total PAR absorbed by the plant canopy per m2 ground over one simulation time step, most of it coming from above Aplantpar with a little coming from below Aplantpar due to reflection. Here, we define for convenience a new variable for the total canopy PAR absorption per ground area:

Aground=Aplantpar+Aplantparμmol/m2 ground

By correcting for the leaf area index Lai (m2 leaf/m2 ground), we get the absorption per leaf area:

Aleaf=AgroundLaiμmol/m2 leaf

However, the leaf photosynthesis model takes incident PAR Fleaf (μmol/m2 leaf/s) as a parameter, not the absorbed PAR. Fleaf is used in the leaf synthesis model to calculate the electron transport rate J (μmol/m2 leaf/s),

J=aFleaf+Jmax(aFleaf+Jmax)24aFleafθJmax2θ

where Jmax (μmol/m2 leaf/s) is the the maximum rate corrected for leaf temperature, and other parameters are as specified for the leaf photosynthesis model. We can derive Fleaf from Aleaf by correcting for leaf PAR absorptivity αleafpar[0;1] and the simulation time step Δt (s),

Fleaf=AleafαleafparΔt

The immediate outputs from the leaf photosynthesis model are all per leaf area, e.g., the net photosynthetic rate Pn (μmol/m2 leaf/s). These rates can be turned into per ground area units by multiplication with the leaf area index to express the rates at canopy level, e.g., Pncanopy=LaiPn (μmol/m2 ground/s).

With a basic parameter setting for the canopy model, Lai=1 and αleafpar=1, the leaf and canopy photosynthesis models give the exact same numerical result, even though one is per leaf area and the other is per ground area. Below (produced by the photosynthesis.R script) the canopy photosynthesis Pncanopy is shown with the same parameter setting for the leaf area model as in the previous figure and with Lai=3 and αleafpar=0.9 in the canopy model:

The self-shading of the canopy brings light intensity down by a factor of Lai which shifts the light response away from asymptotic saturation towards the linear part of the response. This is a simplification, which by itself leads to an overestimate of photosynthesis, since in reality the top leaves might approach light saturation. On the other hand, the plant can compensate for the dwindling light level down through the canopy by producing shade-adapted leaves, i.e. leaves with a light response curve that rises more steeply at low light levels.

Another factor that creates uncertainty in the photosynthesis model is the uncertainty of the boundary layer resistance (rb=gb1). It is difficult to estimate empirically under greenhouse conditions, where it depends on plant architecture as well as air movement caused by ventilation and fans. Thus top leaves are likely to experience a smaller rb than leaves further down in the canopy.

The gradients of light and rb in the plant canopy affects the photosynthesis model in opposite directions. The light gradient can be taken into account by Gaussian integration (Goudriaan, 1986), while the rb gradient must be estimated from published studies. In any case, if the estimation of primary production (e.g. as biomass produced per day) is important, the model will need to be calibrated to the specific greenhouse setting. This includes rb, as well as the species-specific parameters of the leaf photosynthesis model.

If the model is used to explore photosynthesis under different scenarios (e.g., different climate control strategies), the estimated photosynthesis can be used to compare their relative performance. However, the uncertainties and unknowns should be taken into account. The installation of fans will reduce rb, but by how much? Greenhouse temperature influences photosynthesis but have we got accurate estimates of the temperature response of Jmax and Vcmax in the leaf photosynthesis models? One should interpret model predictions of photosynthesis with care.

Climate control

Traditionally, greenhouse climate is controlled by two temperature setpoints: one for ventilation Svent (°C) and one for heating Sheat (°C):

The humidity is controlled indirectly through manipulation of the temperature setpoints, e.g., if it is too humid then the ventilation setpoint is lowered.

In the model, greenhouse air temperature Tin (°C) is considered to be within bounds if it is between the two setpoints ± a tolerance ϵtol (°C). Otherwise, it is either too hot or too cold. The model emulates the climate control computer in a real greenhouse through the logic below, which has five possible outcomes: Turn heating up or down, turn ventilation up or down, or keep current heating and ventilation.

The default corrections of pipe temperature ΔTpipe and ventilation rate Δv should accomodate to the simulation time step Δt (s) :

ΔTpipe=T˙pipeΔtΔv=v˙Δt

where the default rates of change are T˙pipe (°C/s) and v˙ (h1s1).